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Abstract 

In  this  report  a  fast  optimal  nonlinear  filter  is  developed.  This  fast  filter  is  based  on 
two  techniques:  (i)  splitting  of  the  convection  and  diffusion  operators  and  (ii)  tracking  of 
the  “important”  domains  (windows)  with  iterative  reduction  of  the  window  size  -  a  domain 
pursuit  technique.  As  a  result,  the  domain  of  interest  is  determined  adaptively. 

The  developed  nonlinear  filtering  algorithm  is  then  applied  to  a  real  problem  of  tracking 
ballistic  targets  in  six  dimensions.  Simulation  results  for  this  problem  demonstrate  the  fairly 
high  statistical  accuracy,  efficiency,  and  real-time  performance  of  the  proposed  algorithm. 
These  results  also  show  that  the  developed  method  is  much  more  accurate  compared  to  the 
traditional  extended  Kalman  filter. 


1  Introduction 

In  applications  the  target  tracking  problem  can  be  naturally  formulated  as  a  filtering  problem  for 
hidden  Markov  models:  given  an  unobserved  (hidden)  state  (or  signal  or  system)  process,  which 
is  usually  assumed  Markovian,  and  an  observation  (or  measurement)  process,  which  provides 
noisy  information  about  the  state  process,  one  needs  to  estimate  the  state  or  a  function  of  the 
state  at  a  given  time  moment  by  using  all  the  observational  information  available  up  to  that  time 
moment. 

Ideally,  the  involved  stochastic  dynamics  are  linear  and  Gaussian  and  in  this  case  the  filtering 
problem  is  solved  by  the  Kalman  and  Kalman-Bucy  filters.  For  linear  and  Gaussian  models  the 
Kalman  filter  is  optimal  in  the  mean-square  sense  and  had  a  big  success  in  a  wide  variety  of 
applications.  However,  many  real-world  problems  do  not  fit  well  with  linear  dynamic  models. 
Sometimes  one  can  explicitly  describe  the  distribution  of  the  state  given  measurements  (pos¬ 
terior  distribution)  but,  outside  the  realm  of  the  linear  theory,  only  a  very  few  examples  have 
explicitly  described  posterior  distributions.  Since  most  real  problems  are  nonlinear,  this  creates 
a  fundamental  problem.  Successive  linearization  in  short  time  intervals,  the  Extended  Kalman 
Filtering  procedure,  may  be  applied  but  its  serious  disadvantage  is  that  it  often  gives  erroneous 
answers  and  refining  of  computational  effort  can  increase  them. 

Theoretical  study  of  the  general  nonlinear  filtering  problem  has  also  gone  through  more  than 
three  decades  of  efforts  of  mathematicians,  statisticians,  and  engineers  and  has  now  become  more 
or  less  mature  as  a  research  field.  See  the  books  by  Stratonovich  [28],  Jazwinski  [11],  Liptser  and 
Shiryayev  [16],  Kallianpur  [13],  Rozovskii  [26],  Pardoux  [22],  Bensoussan  [2],  Tanizaki  [29],  Elliott, 
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Aggoun,  and  Moore  [4],  In  contrast  to  the  linear  case  where  there  exists  a  finite  dimensional 
statistic,  a  general  nonlinear  problem  is  infinite-dimensional  in  nature.  This  was  the  main  reason 
why  optimal  nonlinear  filters  have  not  been  widely  used  in  real  applications. 

Approximations  to  the  optimal  nonlinear  filter  must  be  adopted.  The  simplest  and  most  widely 
used  approximation  is  the  extended  Kalman  filter  (EKF),  which  is  basically  the  Kalman-Bucy 
filter  applied  to  a  dynamically  linearized  system.  The  EKF  has  been  modified  for  different 
purposes  and  has  many  versions  such  as  the  second  order  approximation  [21]  and  the  iterative 
extended  Kalman  filter  [1,  11,  12]. 

In  addition  to  the  EKF  and  its  modifications,  there  are  two  major  approaches  to  nonlinear  filter 
approximation.  One  approach  assumes  that  the  filtering  densities  belong  to  a  certain  class  of 
functions  such  as  Gaussian  or  exponential  or  some  combinations.  Actually  EKF  is  an  example  of 
an  assumed-density  filter  and  is  perhaps  the  simplest.  Other  examples  are  the  Gaussian  mixture 
filter  (see  [27])  and  more  recently  the  projection  filter  (see  [3]). 

The  other  approach  is  to  use  a  direct  (analytical  or  numerical)  approximation  to  the  optimal  non¬ 
linear  filter.  One  recent  advance  in  this  direction  is  the  Wiener  chaos  decomposition  or  spectral 
separation  scheme  (S'3)  for  nonlinear  filtering  in  continuous  time  [17,  18]  and  similar  (related) 
algorithms  that  also  use  the  “off-line/on-line  separation”  idea  [19,  20].  A  direct  numerical  ap¬ 
proximation  to  the  optimal  nonlinear  filter  is  based  on  computing  the  convolution  integral  in 
the  discrete  filtering  model  [15,  29],  on  using  fast  solvers  for  the  Fokker-Planck  equation  in  the 
continuous-discrete  filtering  model  [14,  19,  24],  or  on  solving  the  Zakai  equation  in  the  case  of 
continuous  time  [6,  8,  9,  10]. 

Both  approaches  encounter  computational  difficulties  in  practical  applications.  The  problem  with 
EKF  and  its  modifications  or,  in  general,  with  the  assumed  density  filters,  is  that  they  do  not 
work  well  if  the  posterior  distribution  differs  from  the  assumed  form.  For  example,  if  one  assumes 
a  Gaussian  distribution  and  uses  EKF  while  in  reality  it  is  multi-peak  (far  from  Gaussian),  then 
the  filter  completely  fails.  Assumed  density  filters  (including  EKF)  fail,  for  example,  in  many 
important  situations  such  as  angle-only  target-tracking  due  to  divergence,  instability,  inaccuracy, 
etc.  The  reason  is  that  the  prespecified  density  class  is  too  restrictive  in  the  general  nonlinear 
case.  Direct  approximation  is  much  better  in  this  class  of  situations  but  has  another  important 
limitation  -  the  “curse  of  dimensionality” .  If,  for  instance,  we  have  a  six-dimensional  model, 
which  is  typical  for  radar  applications,  and  100  points  are  used  in  each  component  of  the  state 
(in  many  cases  100  points  could  even  be  too  few  for  a  satisfactory  estimate  of  the  state),  then 
the  total  number  of  spatial  points  is  N  =  1012.  If  one  uses  a  fast  solver  with  FFT  which  has 
complexity  CdN (\og2  N)d~l  for  dimension  d,  then  one  needs  to  perform  1020C6  flops  at  each 
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time  step,  which  is  simply  unacceptable.  Even  if  the  Fokker-Planck  solver  has  optimal  (linear) 
complexity  (like  ADI),  we  still  need  to  perform  1012C6  flops  each  step.  Note  that  the  constant 
Cd  increases  with  dimension  d  and  C\  >  8. 

What  do  we  do  then?  From  the  one  hand  to  avoid  the  difficulties  that  EKF  faces,  it  is  desirable  to 
use  a  direct  approximation  (rather  than  an  assumed  density  approach)  to  the  optimal  estimate. 
From  the  other  hand  in  order  to  implement  direct  approximations  to  high-dimensional  real 
problems,  we  need  (1)  to  develop  fast  Fokker-Planck  equation  solvers  with  linear  complexity, 
and  (2)  to  reduce  the  number  N  of  spatial  points  substantially  without  reducing  the  accuracy. 
The  former  can  be  achieved  by  the  so-called  operator-splitting  method.  The  latter  goal  is  achieved 
by  reducing  adaptively  at  each  time  step  the  size  of  the  spatial  domain  in  which  the  equation  is 
solved  (the  domain  should  be  reduced  in  each  direction).  This  idea  leads  to  the  adaptive  domain 
pursuit  technique  (DPT)  that  is  developed  below. 

It  is  worth  mentioning  that  the  proposed  domain  pursuit  approach  has  some  similarity  with 
the  EKF.  In  fact,  at  each  time  step  the  DPT  tracks  a  domain  (window)1  in  which  the  target 
is  located  and  then  proceeds  the  nonlinear  filtering  in  this  moving  window  (or  multi- windows), 
whereas  the  EKF  tracks  only  two  parameters  at  each  time  step  (the  mean  and  variance  of  the 
target  state)  and  then  linearizes  the  nonlinear  dynamics  around  the  estimated  mean. 

The  remainder  of  the  report  is  organized  as  follows.  In  Section  2  we  outline  the  basic  facts 
from  the  theory  of  optimal  nonlinear  stochastic  filtering  for  continuous-discrete  time  model. 
Section  3  is  devoted  to  the  development  of  the  fast  algorithm  which  is  based  on  the  convection- 
diffusion  splitting  and  the  domain  pursuit  technique.  In  Section  4  we  apply  the  developed 
general  algorithm  to  a  practically  important  radar  target  tracking  problem  in  six  dimensions. 
The  numerical  results  obtained  in  computational  experiments  are  given  in  Section  4.2.  These 
results  demonstrate  fairly  high  accuracy  and  efficiency  of  the  method  and  show  that  it  is  much 
more  accurate  compared  to  the  EKF. 


2  Continuous-Discrete  Filtering 

2.1  Statement  of  the  problem 

We  are  interested  in  the  continuous-discrete  filtering  model,  since  this  is  perhaps  the  most  ap¬ 
propriate  model  in  real  target  tracking  problems.  Specifically,  the  dynamics  of  target  trajectories 


1This  domain  may  be  multiply  connected  -  multi- windows. 
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is  naturally  continuous  while  the  observations  are  usually  taken  at  discrete  time  moments. 
Consider  the  dynamic  system  described  by  the  stochastic  differential  equation 

dXt  =  b(Xt)dt  +  o(Xt)dWt,  t>  0, 

Xo  ~  7Tq  . 

The  discrete-time  noisy  observations  are  given  by 

Yk  =  hk(Xtk)  +  Rk(Xtk)Vk,  k  =  0,1,...,  (2) 

where  b  :  ]Rd  — >  JRd  is  a  vector-valued  function,  o  :  lRd  — >■  JRdxd  is  a  matrix  with  function  entries, 
{Wt}t>o  is  a  standard  d-dimensional  Brownian  motion  (Wiener  process),  4  =  kr  (r  >  0),  n0  is  a 
(prior)  distribution  of  the  initial  condition,  hk,  Rk  are  given  functions,  and  14  are  i.i.d.  random 
variables.  Without  loss  of  generality,  X0,  {Wt}  and  {14}  are  assumed  to  be  independent,  and 
the  following  regularity  conditions  on  the  parameters  of  the  model  are  also  assumed:  (1)  the 
functions  bk,  hk,  Qk,  Rk,  and  ttq  have  bounded  derivatives  up  to  an  appropriate  order,  and  (2)  all 
the  derivatives  of  7r0  decay  at  infinity  faster  than  any  power  of  \x\. 

For  simplicity,  we  only  consider  the  case  where  the  functions  b  and  o  are  time-independent.  But 
the  discussions  that  follow  can  be  easily  generalized  to  cover  the  time-dependent  case.  In  fact, 
only  the  coefficients  in  the  Fokker-Planck  equation  (3)  below  and  the  associated  semigroup  will 
need  to  be  accordingly  modified. 


2.2  Fokker-Planck  equation 


The  theory  of  stochastic  differential  equations  tells  us  (see  [5,  7,  25])  that  under  certain  conditions 
on  b  and  o,  there  exists  a  unique  solution  Xt  of  (1)  (in  the  sense  of  Ito)  and  that  the  probability 
density  u(t,  x )  of  this  diffusion  process  Xt  satisfies  the  Fokker-Planck  equation  (also  known  as 
the  Kolmogorov  forward  equation) 


du(t,  x ) 
dt 


1  d 
=  E 


2  „™i 


q2  d  q 

(aw(a;)u(i,a;))  -  ^—(bu(x)u(t,x)) 

v=l 


(3) 


where  afiu(x)  is  the  /i-th  row  and  z'-th  column  entry  of  the  product  matrix  cr(x)a(x)*,  and  b„(x) 
is  the  zv-th  component  of  the  vector  b(x). 


Let  T(t)  denote  the  semigroup  associated  with  the  above  Fokker-Planck  equation.  Then  its 
solution  u(t,x )  with  the  initial  value  u(0,  x)  =  g(x)  is  [T(t)g](x). 
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We  now  note  that  since  there  is  no  observation  available  between  tk~ 1  and  4>  the  (prior)  transition 
density  is  the  only  available  information  on  (4-i>4)  and  so  one  has 

p(X,t  =  x  I  Y1-1)  =  [n r)p(Xh_,  =  ■  I  y‘-')]  (*) . 


2.3  Unnormalized  filtering  densities 

By  using  the  Bayes  formula  and  the  Fokker-Planck  equation,  one  can  obtain  the  conditional 
densities  recursively  as  follows 

p(X,„  =  X  I  Yk)  =  [T(t)p(V,,_,  =  •  I  y^1)]  (*),  k  >  1 

p(x 0  =  x\Y°)  =  -^a0{x)Mx) 

where  c{k )  is  the  normalizing  constant  to  make  the  integral  of  p(Xtk  =  x  \  Yk)  to  be  one,  and 
the  “correction  term”  oik{x )  related  to  the  observation  is  given  by 

a„(x)  =  exp  {-i(n  -  AtWnXtMBtWT'tn  ~  h*(*))}  ■  (4) 

The  calculations  are  simplified  if  in  place  of  the  usual  (normalized)  filtering  density  one  uses  the 
unnormalized  filtering  densities  (UFD).  We  define  the  UFD  by2 

Pk(x)  =  (x) [T (r )p*_i] (x) ,  k>l  ^ 

Po(x)  =  a0(x)7r0(x). 


It  is  a  standard  fact  that  for  any  function  g  such  that  E\g(Xt)\  <  oo,  the  conditional  expectation 
of  g(Xtk)  given  Yk  =  a(Y0,  Yu  •  •  • ,  Yk)  can  be  obtained  by 


E(g(Xtk)  |  Yk)  = 


/  .  g(x)pk{x)dx 
J]Rd _ 

J^dpk{x)dx 


(6) 


This  conditional  expectation  is  the  best  mean-square  estimate  of  g(Xtk)  if  E\g(Xt J|2  <  oo. 
(This  is  the  case  if  g  satisfies  |p(x)|  <  K(  1  +  |m|A),  Va;  6  JRd,  for  some  A  and  K  >  0;  see  [16].) 

2 When  Rk  does  not  depend  on  x,  our  definition  here  is  slightly  different  from  the  usually  defined  UFD.  The 
difference  is  in  ak(x),  where  we  keep  the  term  -Yk(RkR*k)~1Yk/2.  This  is  done  to  avoid  computational  instability: 
from  (4)-(5)  we  have  ||pfc||oo  <  ||T(r)p*;_1||00,  Vfc.  In  fact,  the  Fokker-Planck  equation  for  this  UFD  can  also  be 
modified  for  the  purpose  of  stability.  See  Section  2.4.3  and  Section  2.5.2  of  [23]  for  details. 
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3  Adaptive  Domain  Pursuit  Method 

In  this  section,  we  present  a  windowing  technique  which  is  based  on  the  framework  of  splitting 
the  convection  (drift)  and  diffusion  (noise)  operators. 


3.1  Splitting  of  convection  and  diffusion  processes 


To  compute  the  unnormalized  filtering  density,  a  fast  Fokker-Planck  solver  is  needed.  Our  method 
is  based  on  the  operator-splitting  technique. 


Assuming  that  in  the  noise  term  of  (1)  the  covariance  matrix  o  is  constant  and  diagonal,  the 
Fokker-Planck  equation  (3)  becomes 


du(t,x)  _  1  A 
dt  ~  2  dxl 


q2  n  q 

1>—1  uXjs  v 


To  proceed  the  splitting  of  convection  and  diffusion  terms,  denote  by  Tc(t)  and  Td(t)  the  solution 
operators  of  the  equations 


dv(t,x) 

dt 


i/=l 


and 


dw(t,x)  1  A  d2  /  ,  \ 


dt  %  TTi  dxv 

respectively.  Then  it  can  be  proved  that  the  following  approximation  formulae  hold  (see  [23]): 


T(nr)^  =  (Td(r)Tc(r))>  +  0(r), 

=  (Tc(^)Td(T)Tc^))\  +  0(r2). 


(7) 

(8) 


Therefore,  instead  of  solving  the  original  Fokker-Planck  equation,  we  only  need  to  solve  two 
simpler  equations,  for  which  methods  with  linear  computational  complexity  exist. 

We  remark  that  a  big  part  of  computation  in  solving  the  two  simpler  equations  can  be  performed 
before  the  observations  become  available.  This  pre-calculation  substantially  speeds  up  the  on-line 
part  of  the  algorithm. 
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3.2  Domain  pursuit  tracker 

Even  though  we  have  developed  optimal  solvers  for  the  Fokker-Planck  equation,  the  nonlinear 
filtering  problem  can  still  hardly  be  solved  in  real  time  (especially  for  large  state  dimensions). 
The  following  natural  step  is  to  narrow  down  the  size  of  the  domain  in  which  the  Fokker-Planck 
equation  is  solved.  This  can  be  achieved  by  a  windowing  technique  based  on  the  convection- 
diffusion  splitting  framework. 

A  numerical  approximation  to  the  optimal  filter  can  be  expressed  in  the  following  form 

pi  =  Mxi(k))  2  PkA- 1>  k  ^ 1  > 

=  Q'O(a;i(0))7ro(a;i(0)) , 

where  {Pk^}  is  a  (possibly  indirect)  approximation  to  the  fundamental  solution  T(r)  at  step  k. 
Note  that  we  use  the  notation  x*(k)  instead  of  xk  for  the  spatial  points  because  these  are  vectors 
in  Rd  and  x\ ,(k)  will  denote  the  i^-th  component. 

We  remark  that  in  real  implementation,  the  approximate  solution  Shck- i  (or  fikPk-i)  of  the 
Fokker-Planck  equation  is  not  necessarily  computed  directly  from  the  above  summation.  For 
example,  if  an  implicit  scheme  is  used,  then  Sh  contains  an  inverse  matrix  which  is  not  inversed 
directly.  Roughly,  if  the  final  Sh  is  sparse,  then  the  summation  can  be  used  directly;  otherwise 
the  summation  will  take  0(N2)  FLOPS  per  time  step  and  so  some  indirect  technique  such  as 
ADI  should  be  used. 

Assume  we  are  given  an  initial  set 

D0  =  {x1(0),---,xN°(0)} 

of  “important”  points.  This  set  is  chosen  according  to  the  initial  filtering  density.  For  example, 
these  points  can  be  related  to  the  largest  values  of  p0(x)  (or  with  the  most  important  information 
on  Po{x)).  Below  we  describe  how  to  efficiently  compute  the  “important”  points  in  all  the 
subsequent  time  steps  and  also  the  corresponding  values  of  the  unnormalized  filtering  densities 
at  those  points. 

If  the  number  N0  is  too  large,  it  can  be  reduced  in  the  first  several  filtering  steps.  Let  K  >  0  be 
a  small  integer  and  {Nk,  k  >  1}  be  a  decreasing  sequence  of  integers  with  Nk  =  NK  for  k  >  K, 
i.e.  No  <  Ni  <  •••  <  Nk  =  Nk+i  =  •  •  •  Let  L  be  a  positive  integer.  We  will  first  construct  an 
enlarged  set  of  LNk- 1  candidates  for  the  important  points  at  each  step  and  then  choose  from 
them  the  “best”  Nk  points  according  to  the  correction  term  ak. 
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For  simplicity,  we  assume  again  that  the  covariance  matrix  o  is  constant  and  diagonal.  Our 
algorithm  can  be  described  as  follows. 

The  DPT  Algorithm 

•  Initialization:  start  with  time  k  =  0,  domain  D0  of  size  N0,  and  pj  =  p0(V(0)),  1  <  i  <  N0. 

•  Iteration:  for  k  <  K,  run  algorithm  DPT(k)  (see  below)  with  reduction  of  domain  size  Nk. 

•  For  k  >  K ,  run  algorithm  DPT(k)  with  moving  domain  Dk  of  fixed  size  Nk  =  Nk- 1- 

The  algorithm  DPT(k)  proceeds  as  follows 

(1)  for  i  =  1,  •  •  • ,  7Vfc_i,  solve  the  stochastic  differential/integral  equation 

Xt  =  xl(k  -  1)  +  f  b(Xs)ds  +  [  crdWs,  t  £  (0,  A t] , 

with  L  different  sample  paths  of  the  Wiener  process  Wt,  including  the  trivial  case  Wt  =  0, 
and  denote  the  solution  XT  at  time  t  -  r  with  the  j-th  sample  path  of  Wt  by  C'j,  j  = 

i)  •  •  •  >£; 

(2)  determine  the  set  Dk  of  Nk  important  points  x1(k),---,  xNk  ( k )  as  those  with  the  largest 
values  of  ak(C,j)  for  all  i  and  j,  or,  for  each  i,  xl(k)  has  the  largest  value  of  otk(C’j)  for 

j  ~  1)  ‘  ‘  ‘ » 

(3)  for  i  =  1,  •  •  • ,  Nk,  compute 

Pi  =  Mxi(k))  Y,  Ptfpi-i  > 

where  is  computed  according  to  one  of  the  following  two  formulae  which  follow  from 
the  operator-splitting  schemes  discussed  in  the  previous  section: 

(xj{k)  -  xj(k  -  1)  -  Kjx^k  -  1  ))t)2 
2  <72ut 

W  =  exp  |-  f  (*I-W  -g*  -  l))2  _  (v  .  J  , 

and  4  =  {j:l<j<  min(/3^J,p(:_1)  >  r},  r  being  a  thresholding  tolerance. 
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We  stress  that  in  general  the  domains  Dk  may  be  multiply  connected,  i.e.  may  contain  multiple 
windows  (pieces).  The  block  diagram  of  the  domain  pursuit  tracker  is  shown  in  Figure  1. 


Figure  1:  Domain  pursuit  tracker 


The  remarkable  feature  of  the  proposed  filtering  algorithm  is  that  the  domain  of  interest  adap¬ 
tively  changes  and,  as  a  result,  the  number  of  spatial  points  in  the  domain  is  usually  reduced  to 
a  relatively  small  number.  Hence  the  computational  complexity  is  reduced  tremendously  when 
compared  with  the  case  where  a  fixed  size  domain  is  used  for  the  whole  computation.  In  general, 
the  number  N  of  spatial  points  is  exponential  in  d,  i.e.,  N  ~  nd,  but  we  are  reducing  this  number 
for  each  n  of  the  d  dimensions.  Therefore,  the  computational  cost  is  exponentially  reduced  in 
our  algorithms  in  high  dimensions. 


4  Application  to  the  Problem  of  Ballistic  Target  Tracking 

4.1  The  tracking  problem 

To  illustrate  the  performance  of  the  proposed  algorithm,  let  us  consider  a  real  RADAR  tracking 
problem,  which  is  well-known  to  be  difficult  despite  it  does  not  contain  perturbations  in  dynamics. 
(In  principle  we  can  also  handle  a  similar  problem  with  infra-red  or  other  kinds  of  angle-only 
measurements  and  the  model  that  would  include  dynamics  noise.) 
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There  is  a  radar  with  geodetic  latitude  9,  longitude  A,  and  height  h,  observing  a  ballistic  missile 
and  generating  radar  range  R,  azimuth  A,  and  elevation  E  every  r  seconds  (see  Figure  2). 


Figure  2:  Radar  geometry 


The  missile  is  assumed  to  be  in  unpowered  ballistic  flight  whose  six  dimensional  dynamic  equa¬ 
tions  of  motion  are 


dX(t) 

dt 


=  V(t) , 


dV(t)  _  X(t) 

dt  ~  Vwir 


t  >  0, 


(9) 


where  p  =  3.986012  x  1014,  X(t)  =  (Xi(t),X2(t),X3(t))  and  V(t)  =  (Vi (t),V2(t),V3(t))  are  the 
position  and  velocity  of  the  missile  at  time  moment  t,  and  X(0)  and  V(0)  are  Gaussian  random 
vectors  with  known  mean 


E{X{ 0))  =  [0,0,7.45005724  x  106]T(m), 

E(V{ 0))  =  [-3.96745  x  103,  -2.37208  x  103, 2.15685  x  103]T (m/sec), 


and  covariance 


Cov(X(0),  17(0))  =  diag[4  x  106, 4  x  106, 4  x  106, 104/3, 104/3, 104/3] . 


The  computation  of  the  radar  measured  range,  azimuth,  and  elevation  from  the  missile’s  true 
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inertial  position  at  time  tk  =  hr  proceeds  as  follows: 

R(k)  =  \\XL(tk)\\  +  £ivi(k) , 

A(k)  =  tan-1  +  £2^2 (&) , 

m=™'wM\+eMk)- 


where  S\  =  16,  £2  =  £3  =  \/3. 04617  x  10  6,  and 


(10) 


XL(t)  =  Tloc,ecbf  •  Pecef,eci(£)  •  X(t)  —  As] . 


Here  LOC,  ECEF,  and  ECI  stand  for  the  three  relevant  coordinate  systems:  ECI  for  Earth 
Centered  Inertial  True-of-Date,  ECEF  for  Earth  Centered  Earth  Fixed,  and  LOC  for  Local 
East-North-Up.  The  coordinate  transformation  matrices  Tloc,ecef  and  2ecef,eci(^)  and  the 
radar  site  location  Xs  are  given  by 


As 


TLoc.ecef  — 


Tecef.eciW 


—  sin  A  —  sin  9  cos  A  cos  9  cos  A 

cos  A  —  sin  9  sin  A  cos  9  sin  A 

0  cos  9  sin  9 

cos  (ujet)  sin(o;et)  0 
=  —  sin(a>et)  cos  (a iet)  0  , 

0  0  1 


5 


(i a  +  h)  cos  9  cos  A 
(a  -1-  h)  cos  9  sin  A  , 
(P  +  h )  sin  9 


Q>e 

\f\  —  e2  sin2  9  ’ 


P  = 


ae(l  -  e2) 

\/l  —  e2  sin2  9  ’ 


where  9  =  1.12032684685  (rad),  A  =  -2.60246044764  (rad),^  =  7.2722052162296xl0~5(rad/sec), 
ae  =  6.37815  x  106  (m),  e2  =  6.69342162296  x  10~3,  h  =  0  (m). 


4.2  Monte  Carlo  simulations 

Below  we  present  the  results  of  computational  experiments.  (L  =  1  was  taken  for  Algorithm  1.) 

Assume  the  observations  are  available  at  every  r  =  1  second.  Since  there  is  no  noise  in  the  state 
dynamics  (9) ,  we  do  not  need  the  diffusion  smoothing  as  described  for  the  general  model  in  the 
algorithm.  If  the  noise  is  added  in  the  dynamics,  then  the  general  algorithm  should  be  applied. 

In  our  simulation,  we  ran  our  filter  200  times  with  random  initial  conditions  and  random  obser¬ 
vations.  We  used  two  different  values  of  N,  the  number  of  spatial  points  in  the  moving  domain 
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Table  1:  Computational  complexity  of  the  algorithm 


Simulation 

N 

flops 

cpul 

cpu2 

Simulation  1 

729 

72,182 

0.14 

0.04 

Simulation  2 

15,625 

1,546,886 

0.64 

0.48 

(window).  Specifically,  N  =  36  =  729  and  N  =  56  =  15, 625.  The  average  computational  cost 
for  the  two  experiments  in  terms  of  CPU  seconds  and  FLOPS  per  time  step  are  given  in  Table  1, 
where  cpul  and  cpu2  are  the  CPU  time  for  one  step  of  calculations  with  and  without  graphics, 
respectively.  The  computation  was  performed  with  Matlab  on  a  Sun  Ultra  Enterprise  4000  at 
the  University  of  Southern  California.  From  this  table,  it  is  clear  that  real-time  performance  has 
been  achieved  even  for  tracking  10  to  20  targets. 

For  both  values  of  N,  we  tested  three  different  situations.  In  the  first  case  we  assume  that  the 
true  initial  location  of  the  target  and  its  velocity  are  exactly  on  the  grid;  in  the  second  case  it 
is  assumed  only  that  the  initial  velocity  is  exactly  on  the  grid;  and  in  the  third  case  both  the 
initial  location  and  velocity  are  not  on  the  grid  (the  most  realistic  case).  Of  course,  in  either 
case,  the  filter  does  not  know  the  true  initial  location  of  the  target.  In  the  first  case  performance 
is  perfect:  with  N  -  729,  the  average  errors  for  both  X (t)  and  V (t)  became  close  to  zero  after 
t  =  140  seconds.  The  average  errors,  defined  as 

Xerr,«  =  ^E  |  and  7„,(i|  =  —  £  |  V^(t)  -  KW«)I 

( y  =  1,2,3)  are  shown  in  Figure  4  (for  the  two  simulations).  Note  that  because  the  initial  error 
in  X(0)  is  too  large,  it  is  not  shown  in  the  picture.  See  Figure  3  for  the  average  errors  (including 
the  initial  errors)  according  to  the  state  dynamics  without  any  observed  information. 

In  the  second  and  third  cases,  when  the  initial  variance  (for  the  true  initial  state)  is  relatively 
small,  we  also  obtained  good  results.  For  example,  in  the  second  case,  with  the  initial  variance 

Cov(X(0),  U(0))  =  diag[104, 104, 104, 10, 10, 10] , 

the  average  errors  in  V  go  to  zero  and  the  average  errors  in  X  begin  to  decrease  after  t  =  90 
seconds.  These  errors  are  shown  in  Figure  5.  Again  the  initial  error  in  X(0)  is  not  shown  in  the 
picture.  In  the  third  case,  with  the  initial  variance 

Cov(X(0),  U(0))  =  diag[104, 104, 104, 1, 1, 1] , 

the  average  errors  behave  similarly  except  that  they  decrease  at  a  slower  rate. 
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The  real  initial  errors  used  in  experiments  are  listed  in  Table  2.  Typical  errors  at  t-  100  seconds 
(with  the  original  large  error)  are  shown  in  Table  3.  This  table  contains  data  for  DPT  and  EKF. 
It  may  be  seen  that  the  proposed  method  is  much  more  accurate. 


Table  2:  Real  initial  errors 


Error  in  X3  (m)  Error  in  X2  <m)  Error  in  XI  (m)  Error  in  X3  (m)  Error  in  X2  (m)  Error  in  XI  (m) 
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M.C.  Errors  in  X  M.C.  Errors  in  V 


Time  t  (sec) 


Figure  5:  Errors  in  Case  2,  N  =  36 


5  Conclusion 

1.  In  this  report  we  described  the  developed  nonlinear  filtering  algorithm  that  is  based  on  the 
“domain  pursuit  method” .  This  method  is  directed  towards  obtaining  robust  nonlinear  tracking 
algorithms  with  manageable  complexity  and  high  statistical  performance  (close  to  the  optimal 
level) . 

2.  The  algorithm  is  applied  to  a  realistic  problem  that  is  typical  for  tracking  ballistic  missiles  by 
radar.  The  considered  scenario  includes  targets  with  “hard”  trajectories  that  should  be  localized 
in  100-150  seconds.  The  results  of  simulation  show  that  the  developed  algorithm  substantially 
outperforms  the  conventional  EKF  tracker  it  terms  of  mean-squared  tracking  error  and  at  the 
same  time  has  satisfactory  computational  complexity  (may  be  applied  in  real  time). 
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